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Abstract 

In this paper we propose a novel two-step linear optimization model to calculate energy-efficient 
timetables in metro railway networks. The resultant timetable minimizes the total energy consumed by 
all trains and maximizes the utilization of regenerative energy produced by braking trains, subject to 
the constraints in the railway network. In contrast to other existing models, which are MV-haid, our 
model is computationally the most tractable one being a linear program. We apply our optimization 
model to different instances of service PES2-SFM2 of line 8 of Shanghai Metro network spanning a full 
service period of one day (18 hours) with thousands of active trains. For every instance, our model 
finds an optimal timetable very quickly (largest runtime being less than 13s) with significant reduction 
in effective energy consumption (the worst case being 19.27%). Code based on the model has been 
integrated with Thales Timetable Compiler - the industrial timetable compiler of Thales Inc that has 
the largest installed base of communication-based train control systems worldwide. 

Keywords Railway networks, energy efficiency, regenerative braking, train scheduling, linear program¬ 
ming. 


1 Introduction 

1.1 Background and motivation 

Efficient energy management of electric vehicles using mathematical optimization has gained a lot of attention 
in recent years HI mu [5]. When a train makes a trip from an origin platform to a destination platform, 
its optimal speed profile consists of four phases: 1) maximum acceleration, 2) speed hold, 3) coast and 4) 
maximum brake [HI as shown in Figure in a qualitative manner. Most of the energy required by the train is 
consumed during the accelerating phase. During the speed holding phase the energy consumption is negligible 
compared to accelerating phase, and during the coasting phase there is no need for energy. When the train 
brakes, it produces energy by using a regenerative braking mechanism. This energy is called regenerative 
braking energy. Calculating energy-efficient timetables for trains in railway networks is a relevant problem 
in this regard. Electricity is the main source of energy for trains in most modern railway networks; in such 
networks, a train is equipped with a regenerative braking mechanism that allows it to produce electrical 
energy during its braking phase. In this paper, we formulate a two-step linear optimization model to obtain 
an energy-efficient timetable for a metro railway network. The timetable schedules the arrival time and the 
departure time of each train to and from the platforms it visits such that the total electrical energy consumed 
is minimized and the utilization of produced regenerative energy is maximized. 
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1.2 Related work 


The general timetabling problem in a metro railway network has been studied extensively over the past 
three decades [7]. However, very few results exist that can calculate energy-efficient timetables. Now we 
discuss the related research. We classify the related work as follows. The first two papers are mixed integer 
programming model, the next three are models based on meta-heuristics and the last one is an analytical 
study. 

A Mixed Integer Programming (MIP) model, applicable only to single train-lines, is proposed by Pea- 
Alcaraz et al. [8] to maximize the total duration of all possible synchronization processes between all possible 
train pairs. The model is then applied successfully to line three of the Madrid underground system. However, 
the model can have some drawbacks. First, considering all train pairs in the objective will result in a 
computationally intractable problem even for a moderate sized railway network. Second, for a train pair in 
which the associated trains are far apart from each other, most, if not all, of the regenerative energy will be 
lost due to the transmission loss of the overhead contact line. Finally, the model assumes that the durations 
of braking and accelerating phases stay the same with varying trip times, which is not the case in reality. 

The work in proposes a more tractable MIP model, applicable to any railway network, by considering 
only train pairs suitable for regenerative energy transfer. The optimization model is applied numerically to 
the Dockland Light Railway and shows a significant increase in the total duration of the synchronization 
process. Although such increase, in principle, may increase the total savings in regenerative energy, the 
actual energy saving is not directly addressed. Similar to [5], this model too, assumes that even if the trip 
time changes, the duration of the associated braking and accelerating stay the same. 

Other relevant works implement meta-heuristics. The work in m implements genetic algorithm to 
calculate timetables that maximize the utilization of regenerative energy while minimizing the tractive energy 
of the trains. Numerical studies for the model is implemented to Beijing Metro Yizhuang Line of China 
showing notable increase in energy efficiency. The work in m presents a cooperative integer programming 
model to utilize the use of regenerative energy of trains and proposes genetic algorithm to solve it. Similar 
to m, this numerical studies have been performed to Beijing Metro Yizhuang Line of China, though the 
improvement is stated in the increase in overlapping time only. The work in |12| presents a nonlinear 
integer programming model which is solved using simulated annealing. The numerical experiments have 
been conducted for the island line of the mass transit system in Hong Kong. 

An insightful analytical study of a periodic railway schedule appears in m- The model uses the KKT 
conditions to calculate and analyze the properties of an energy-efficient timetable. The resultant analytical 
model is then applied to Beijing Metro Yizhuang Line of China numerically, which shows that the model 
can reduce the net energy consumption considerably. 


First optimization model 


Second ojrtimization model 



Figure 2: Flow-chart of the two steps of the optimization model 
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1.3 Contributions 


We propose a novel two-step linear optimization model to calculate an energy-efficient railway timetable. 
The first optimization model minimizes the total energy consumed by all trains subject to the constraints 
present in the railway network. The problem can be formulated as a linear program, with the optimal value 
attained by an integral vector. The second optimization model uses the optimal trip time from the first 
optimization model and maximizes the transfer of regenerative braking energy between suitable train pairs. 
Both the steps of our optimization model are linear programs, whereas the optimization models in related 
works are JVV-hard. A flow-chart of the two steps of the optimization model is shown in Figurej^ Our model 
can calculate energy-efficient railway timetables for large scale networks in a short CPU time. Code based on 
the model has been embedded with the railway timetable compiler (Thales Timetable Compiler) of Thales 
Inc, which has the largest installed base of communication-based train control systems worldwide. Thales 
Timetable Compiler is used by many railway management systems worldwide including: Docklands Light 
Railway in London, UK., the West Rail Line and Ma On Shan Line in Hong Kong, the Red Line and Green 
Line in Dubai, the Kelana Jaya Line in Kuala Lumpur. 

This paper is organized as follows. In Section we describe the notation used, and then in Section we 
model and justify the constraints present in the railway network. The first optimization model is presented in 
Section]^ Section [^formulates the second optimization problem that additionally maximizes the utilization 
of regenerative braking energy. Section describes the limitations of the optimization model. In Section 
we apply our model to different instances of an existing metro railway network spanning a full working day 
and describe the results. Section presents the conclusion. 


2 Notation and notions 

Every set described in this paper is strictly ordered and finite unless otherwise specified. The set-cardinality 
(number of elements of the set) and the fth element of such a set C is denoted by IC] and C{i) respectively. 
The set of real numbers and integers are expressed by R and Z respectively; subscripts -I- and -|--|- attached 
with either set denote non-negativity and positivity of the elements respectively. A column vector with all 
components one is denoted by 1. The symbol ^ stands for componentwise inequality between two vectors and 
the symbol A stands for conjunction. The number of nonzero components of a vector x is called cardinality 
of that vector and is denoted by card(x). Note that, cardinality of a vector is different from set-cardinality. 
The zth unit vector is the vector with all components zero except for the fth component which is one. 
The epigraph of a function / : C —>■ R (where C is any set) denoted by epi / is the set of input-output pairs 
that / can achieve along with anything above, i.e., 

epi f = {{x,t) G C X 'R \xGC,t> fix)}. 

The convex hull of any set C, denoted by conv C, is the set containing all convex combinations of points in 
C. Consequently, if C is nonconvex, then its best convex outer approximation is conv C, as it is the smallest 
set containing C. 

The set of all platforms in a railway network is indicated by Af. A directed arc between two distinct and 
non-opposite platforms is called a track. The set of all tracks is represented by A. The directed graph of the 
railway network is expressed by (A/”, A). The set of all trains is denoted by T, where |T| is fixed. The sets of 
all platforms and all tracks visited by a train t in chronological order are denoted by Af* C Af and A* A 
respectively. The decision variables are the train arrival and departure times, to and from the associated 
platforms, respectively. Let a* and d\ be the arrival time and the departure time of the train t G T to and 
from the platform i G Af*. 

Table [T] contains the list of symbols used in this paper. 

Table 1: List of symbols 

\C\ The number of elements of a finite countable set C 
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The ith element of a finite countable set C 

The set of real numbers 

The ith unit vector 

The epigraph of a function / 

The number of nonzero elements in a vector x G R” 

The set of all platforms in a railway network 
The set of all tracks 
The set of all trains 

The set of all platforms visited by a train t in chronological order 

The set of all tracks visited by a train t in chronological order 

The arrival time of train t at platform i 

The departure time of train t from platform i 

The trip time window for train t from platform i to platform j 

The set of all train pairs involved in turn-around events on crossing-over {i,j) 

The trip time window for train t on the crossing-over {i,j) 

The dwell time window for train t at platform i 

The set of all platform pairs situated at the same interchange stations 
The set of connecting train pairs for a platform pair (i,j) € x 

The connection window between train t at platform i and and train t' at platform j 
The time bound for any event time in the timetable 
The set of train-pairs who move along that track (i,j) 

The headway time window between train t and t' at or from platform i 

The total travel time window for train t to traverse its train path 

The set of all nodes in the constraint graph 

The set of arcs in the constraint graph 

The set of all arcs associated with trip time constraints 

The arrival or departure time of some train from a platform in the constraint graph 

The time window associated with arc {i,j) of the constraint graph 

Energy consumption associated with the trip (i,j) € .Atrip 

Affine approximation for fij 

Solution to step one optimization model 

The set containing all opposite platform pairs powered by the same electrical substations 
The set of all trains which arrive at, dwell and then depart from platform i 

Temporally closest train to the right of train t 

Temporally closest train to the left of train t 
Temporally closest train to train t 

The relative distance of a* from the regenerative alignment point 
The relative distance of the consumptive alignment point from d* 

The set of all synchronization processes between suitable train pairs 

A subset of £ containing elements of the form t ) 

A subset of £ containing elements of the form t ) 

Convex envelope of function / 


3 Modelling the constraint set 

In this section we describe the constraint set for our optimization model. This comprises the feasibility 
constraints for a railway network of arbitrary topology, and the domain of the decision variables. In every 
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active railway network, the railway management has an operating feasible timetable; we use the sequence of 
the trains from that timetable. The lower and upper bound of the constraints are integers representing time 
in seconds. 

3.1 Trip time constraint 

The trip time constraints play the most important role in train energy consumption and regenerative energy 
production. These can be of two types as follows. 

3.1.1 Trip time constraint associated with a track 

Consider the trip of any train t € T from platform i to platform j along the track (qj) S . The train t 
departs from platform i at time d\, arrives at platform j at time a*, and it can have a trip time between tC 
and rb. The trip time constraint can be written as follows: 

VteT, y{i,j)eA\ < a] - d\<T%. (1) 

3.1.2 Trip time constraint associated with a crossing-over 

A crossing-over is a special type of directed arc that connects two train-lines, where a train-line is a directed 
path with the set of nodes representing non-opposite platforms and the set of arcs representing non-opposite 
tracks. If after arriving at the terminal platform of a train-line, a train turns around by traversing the 
crossing-over and starts travelling through another train-line, then the same physical train is treated and 
labelled functionally as two different trains by the railway management [H page 41]. Let ip be the set of all 
crossing-overs, where turn-around events occur. Consider any crossing-over (i,j) S p, where the platforms 
i and j are situated on different train-lines. Let Bij be the set of all train pairs involved in corresponding 
turn-around events on the crossing-over (z,j). Let S Bij. Train t £ T turns around at platform i 

by travelling through the crossing-over {i,j), and beginning from platform j starts traversing a different 
train-line as train t' £T \ {t}. A time window [k-* ,7t-* ] has to be maintained between the departure of the 
train from platform i (labelled as train t) and arrival at platform j (labelled as train t'). We can write this 
constraint as follows: 


'd{i,j)ep, e B^j, k^'< a*- - dl < . (2) 

To clearly illustrate the constraint we consider Figure]^ Here we have two train lines: line 1 and line 2. The 
terminal platform on line 1 is platform i and the first platform on line 2 is platform j. The crossing-over 
from line 1 to line 2 is the arc {i,j). The train shown in the figure is labelled as t on platform i and labelled 
as train t' on platform j. 

3.2 Dwell time constraint 

When any train t G T arrives at a platform i S Af*, it dwells there for a certain time interval denoted by 
so that the passengers can get off and get on the train prior to its departure from platform j. The 
dwell time constraint can be written as follows: 

VteT, VzGAA‘, ^,<d\-a\<5\. (3) 

Every train t € T arrives at the first platform A/'*(l) in its train-path either from the depot or by turning 
around from some other line, and departs from the final platform A/'*(|A/'*|) in order to either return to the 
depot or start as a new train on another line by turning around. So, the train t dwells at all platforms in 
A/"*. This is the reason why in Equation the platform index i is varied over all elements of the set A/"*. 
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Terminal platform i on line 1 



Figure 3: Trip time constraint associated with a crossing-over 


3.3 Connection constraint 

In many cases, a single train connection might not exist between the origin and the desired destination of 
a passenger. To circumvent this, connecting trains are often used at interchange stations. Let % C A/” x A/” 
be the set of all platform pairs situated at the same interchange stations, where passengers transfer between 
trains. Let Cij be the set of connecting train pairs for a platform pair (i, j) € y. For a train pair (t, t') € Cij, 
train t is arriving at platform i and train t' € 'T is departing from platform j. A connection time window 
denoted by ,Xij ] is maintained between arrival of t and subsequent departure of t', so that passengers 
can get off from the first train and get on the latter. Let {i,j) £ y. Then the connection constraint can be 
written as: 

V(Li) e X, V(Li') e c^-, x“' < d*' - a* < xlj ■ (4) 


3.4 Headway constraint 


In any railway network, a minimum amount of time between the departures and arrivals of consecutive trains 
on the same track is maintained. This time is called headway time. For maintaining the quality of passenger 
service, many urban railway system keeps an upper bound between the arrivals and departures of successive 
trains on the same track, so that passengers do not have to wait too long before the next train comes. Let 
(L j) £ A be the track between two platforms i and j, and Hij be the set of train-pairs who move along that 

track successively in order of their departures. Consider (t, t') £ Hij, and let [hf , /if ] and [df , df ] be the 
time windows that have to be maintained between the departures and arrivals of the trains t and t' from 
and to the platforms i and j respectively. So, the headway constraint can be written as: 


v(i,j)£A, 


Ik' < 4 -d\< hf' 


/if <a‘'-a‘</if. 


(5) 


Similarly, headway times have to be maintained between two consecutive trains going through a crossing 
over. Consider any crossing over {i,j) £ f and two such trains, which leave the terminal platform of a 
train-line i labelled as ti and t 2 , traverse the crossing over (i, j), and arrive at platform j of some other 
train-line labelled as and The set of all such train quartets (t 2 T 2 )) is represented by Hij. Let 

[/if*^,/if be the headway time window between the departures of trains ti and t 2 from platform i and 

[/if /if be the headway time window between the arrivals of the trains t\ and t '2 to the platforms j. The 
associated headway constraints can be written as: 


V(L j) £ V((ii,/'i), (t2T2)) G idt. 




t\t2 


< d? - d7 < K 


-tit2 


A 


-3 


^ ^2 ^ 7^1*2 

< a/ - a/ < 


( 6 ) 
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Now we discuss the relation between passenger demand, headway and number of trains. We denote the 
passenger demand by D, train capacity by c and utilization rate by u. If we denote the number of trains in 
service per hour by n, then we have 

D = c X u X n 

[TU] , Because the headway time h satisfies the relation h = we have 

, 3600 X cx u 


It should be noted that the train capacity c and the utilization rate u are constant parameters. However the 
passenger demand varies with time. As a result, in the equation above trains will have different headway at 
different periods. 

3.5 Total travel time constraint 

The train-path of a train is the directed path containing all platforms and tracks visited by it in chronological 
order. To maintain the quality of service in the railway network, for every train t G 'T, the total travel time 
to traverse its train-path has to stay within a time window [Tp,Tp]. We can write this constraint as follows: 

Vt G T, Tp < < Tp, (8) 

where and are the first and last platform in the train-path of t. 

3.6 Domain of the event times 

Without any loss of generality, we set the time of the first event of the railway service period, which corre¬ 
sponds to the departure of the first train of the day from some platform, to start at zero second. By setting 
all trip times and dwell times to their maximum possible values we can obtain an upper bound for the final 
event of the railway service period, which is the arrival of the last train of the day at some platform, denoted 
by m G Z++- So the domain of the decision variables can be expressed by the following equation: 

Vt G T, Vf G A/”*, 0 < a- < TO, 0 < d- < m. (9) 

In vector notation the decision variables are denoted by a = and d = 


4 First optimization model 

In this section we formulate the first optimization model that minimizes the total energy consumed by all 
trains in the railway network. The organization of this section is as follows. First, in order to keep the proofs 
less cluttered, we introduce an equivalent constraint graph notation. Then we formulate and justify the first 
optimization problem. Finally we show that the nonlinear objective of the initial optimization model can be 
approximated as a linear one by applying least-squares. This results in a linear optimization problem, which 
has the interesting property that its optimal solution is attained by an integral vector. 

4.1 Constraint graph notation 

Each of the constraints described by Equations 0® is associated with two event times (either arrival or 
departure time of trains at stations), where one of them precedes another by a time difference dictated by the 
time window of that constraint. This observation helps us to convert our initial notation into an equivalent 
constraint graph notation which we describe as follows. 



Converting the initial notation into an equivalent constraint graph 

• Nodes of the constraint graph: All event times in the original notation are treated as nodes in the 
constraint graph, the set of those nodes is denoted by Af and the value associated with a node i G Af is 
denoted by Xi, which represents the arrival or departure time of some train from a platform. Consider 
any two nodes in the constraint graph; if there exists a constraint between the two in the original 
notation, then in the constraint graph we create a directed arc between them, the start node being the 
first event and the end node being the later one. The set of arcs thus created in the constraint graph 
is denoted by A. Note that there cannot be more than one arc between two nodes in the constraint 
graph. 

• Ares of the constraint graph: With each arc {i,j) G .A we associate a time window [lij,Uij] with their 
values determined from the Equations Q-([^. So, each arc (i,j) G A corresponds to a constraint of 
the form hj < xj — Xi < Uij. The set of all arcs associated with trip time constraints is expressed by 
Atrip C A. 



dl al 


Figure 4: Example of a very simple railway network. 


Example Figure [^represents a very simple network with 7 stations represented by black dots and 2 trains 
represented by the directed rectangles. The pointed edge of the symbol indicates the direction. In this 
network, T = {1,2} and and the stations are enumerated as {1,2, 3,4, 5,6, 7}. There are two train lines 
as shown in the figure. It should be noted that nodes represent stations, not platfrms. Nodes 2 and 3 
represent interchange stations, where there are two platforms for both lines on different levels. So, the event 
times at node 2 and node 3 are associated with different platforms and are differentiated using black and 
magenta colors. The set of tracks visited by train 1 is = {(1, 2), (2, 3), (3,4)}, and the set of tracks 
visited by train 2 is A^ = {(5, 2), (2, 6), (6, 3), (3, 7)}. The set of all tracks is then A = A^ U A?. The 
event times corresponding to train 1 and train 2 are shown in black and magenta colours respectively in 
FigureApplying the conversion process described above we can convert the initial notation in Figure]^ 
to the constraint graph shown in Figure . 
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4.2 Formulation of the first optimization model 

A train consumes most of its required electrical energy during the acceleration phase of making a trip 
from an origin platform to a destination platform. Trip time constraints play the most important role in 
energy consumption and regenerative energy production of trains. Once the trip time for a trip is fixed, 
an energy optimal speed profile can be calculated efficiently by existing software [T5] . [HI page 285], such 
as Thales Train Kinetics, Dynamics and Control (TKDC) Simulator in our case. The TKDC simulator 
assumes maximum accelerating - speed holding - coasting - maximum braking strategy for calculation of 
speed profile. Theoretically this is the optimal speed profile according to the monograph |5]. For calculation 
of the optimal speed prohle of a train while making a trip, we refer the interested reader to the highly cited 
papers na im m m HD]. An excellent state-of-the art review of calculating energy-efficient or energy- 
optimal speed profile can be found in m- The key assumptions that we have made in calculating the speed 
profiles, all of which are well-established in relevant literature |2T], are as follows. 

• The acceleration and braking control forces are assumed to be uniformly bounded or bounded by 
magnitude constraints. These magnitude constraints are assumed to be dependent on the speed of the 
train and non-increasing. 

• The initial speed of the train from the origin platform and final speed at the destination platform are 
assumed to be zero. The speed is assumed to be strictly positive at every other position. 

• Only the force associated with positive acceleration consumes energy. 

• The rate of energy dissipation from frictional resistance is assumed to be a quadratic function of the 
form oo -l- aiv + a 2 V^ (also known as the Davis formula [22] 1. where v stands for speed of the train 
and ao,ai,a 2 are non-negative physical constants. These constants can be determined approximately 
in practice by measuring the difference between the applied tractive force and the net tractive force. 

• We have assumed the position of the train to be an independent variable and speed of the train to be 
a dependent variable. This assumption appears all the major contributions. This assumption appears 
all the major contributions. A notable exception is |19| . which considers kinetic energy and time as 
the dependent variables and position as the independent variable. 

The electrical power consumption and regeneration of a train on a track is determined by its speed profile, 
so the optimal speed profile also gives the power versus time graph {power graph in short) for that trip. 
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However, in the total railway service period there are many active trains, whose movements are coupled by 
the associated constraints. So, finding the energy-minimal trip time for a single trip in an isolated manner 
can result in a infeasible timetable. Consider an arc (*,j) S -dtrip in the constraint graph, associated with 
some trip time constraint. Let us denote the energy consumption function for that trip fij : R+-|- R+-i- 

with argument {xj — Xi). The first optimization problem with the objective to minimize the total energy 
consumption of the trains can be written as: 

minimize “ ^^0 

subject to lij < Xj — Xi < Uij, V(i, j) G A (10) 

0 < Xi < m, yi G Af, 


where the decision vector is {xi)i^f^ G 

The exact analytical form of every component of the objective function, i.e., fij{xj — Xi) for (i,j) G Htop 
is not known and may be intractable [18] . However, irrespective of the exact analytical form, the energy 
function can be shown to be monotonically decreasing in trip time, i.e., it is non-increasing with the increase 
in trip time, if the optimal speed profile is followed [23] . Even when a train is manually driven with 
possibly suboptimal driving strategies, the average energy consumption of the train is found empirically to 
be monotonically decreasing in the trip time [24] . 

Also, the energy function is relatively easy to measure in practice [i Section 1.5]. For any (i,j) € Atrip, 


we denote the measured trip times ..., — xf’) and the corresponding energy consumption 

data 

In any subway system, the amount by which the trip time is allowed to vary in Equations Q and ([^ is 
typically on the order of seconds [25] . which motivates us to make the following assumption. 


yph 


Assumption 1. The amount by which the trip time is allowed to vary is on the order of seconds, i.e., for any 
the trip time window is on the order of seconds. 

The monotonically decreasing nature of the energy function together with Assumption allows us to 
approximate the energy function fij{xj — Xi) as an affine function. Recall that in practice, we can measure 
the energy ..., and associated trip times ..., — x-^^), which is obtainable easily 

with present technology [S] Section 1.5]. 

Now we want to formulate an optimization problem which will provide us with the best possible affine 
approximation of the energy function fij(xj — Xi). We do so by applying least-squares and fit a straight 
line through measured energy versus trip time data. We seek an affine function Cij{xj — Xi) -G bij = {xj — 
Xi, l)'^{cij, bij) where we want to determine Cy and bij. 

The affine function approximates the measured energy in the least-squares sense as follows: 


{c,j,bij) = argmin^ (cij(xf ^ - xf -f kj - 




= argmin 

{cijMj) 




(xj^^ - x^l\ 1)^ 





>(!)■ 

Jlj 



Cij 





pi]. 


1 

_1 



( 11 ) 


The problem above is an unconstrained optimization problem with convex quadratic differentiable objective. 
So it can be solved by taking the gradient with respect to {cij,bij), setting the result equal to zero vector 
and then solving for (c^, bij ). This yields the following closed form solution: 


Cij 


f 

{xf'^ - xf\l)'^ 

T 

(xf^ -xf\l)^ 

\ 

pij 


[ 

(x^-^^ - x-^\ 1)^ 


(xj^^ - x-^\ 1)^ 

/ 


(xf^ -xf\l)^ 


l(Xj-^^ - x[^\ l)"^ 


fW 

Jlj 


f{p) 

^^3 J 


( 12 ) 
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Using Equation (|12p 
an affine one: 


we can approximate the nonlinear objective of the optimization problem (10) as 


j)&Atrip ~ ^j) + ^ measurement of the quality of such httings is given by the 

coefficient of determination, which can vary between 0 to 1, with 0 being the worst and 1 being the best |261 
page 518]. In our numerical studies the mean coefficient of determination of the energy fittings over all the 
different trips of all the trains is found to be 0.9483 with a standard deviation of 0.05, which justifies our 
approach. We can also discard the bijS from the objective, as it has no impact on the minimizer. Thus we 
arrive at the following linear optimization problem to minimize the total energy consumption of the trains: 


minimize 

subject to lij < Xj — Xi < Uij, '^{i,j) G A (13) 

0 < Xi < m, Vi € Af. 


Note that, we have not used the same cost-time curve for all the trips. Each of the constituent parts 
Cij{xj — Xi) of the objective function j)^AtTi ~ Problem (12) represents approximated affine 

function for each of the trips considered in the optimization problem. If optimal speed profiles for the trips 
are available to the railway management, from the first optimization model the optimal trip times for those 
optimal speed profiles can be found. If available speed profiles are suboptimal, then the first optimization 
model would still produce an energy-efficient timetable with best trip times subject to the available speed 
profiles. 

An important property of this optimization model is that the polyhedron associated with optimization 
problem has only integer vertices, so the optimal value is attained by an integral vector. A necessary and 
sufficient condition of integrality of the vertices of a polyhedron is given by the following theorem HZl page 
269, Theorem 19.3], which we will use to prove the subsequent proposition. 


Theorem 1. Let A be a matrix with entries 0,4-1, or —1. For all integral vectors a,b,c,d the polyhedron 
{x G R’" j c ^ a; ^ d, a ^ Ax ^ b} has only integral vertices if and only if for each nonempty collection of 
columns of A, denoted by C, there exist two subsets , Cf and C 2 such that CiU C 2 = C, Ci 0 02 = 0 , and 
the sum of the columns in Ci minus the sum of the columns in C 2 is a vector with entries 0,1 and —1. 


Proposition 1. The optimization problem (13) has an integral optimal solution. 


Proof. We write the problem (13) in vector form. We construct a cost vector c, such that a component of 
that vector is if it is associated with a trip time constraint in the original notation, and zero otherwise. 
Construct integral vectors I = {hj)(ij)^A^ ^ matrix A G {—1,0, such that the 

(k,i)th entry of the matrix A, denoted by aki, is associated with the kth hyperarc and ith node of the 
constraint graph as follows: 


{ 1 if node i is the end node of hyperarc k, 

— 1 if node i is the start node of hyperarc k, 
0 otherwise. 


So, the vector form of the optimization problem (13) is: 


minimize cFx 
subject to I ^ Ax ^ u, 
0 ^ a: ^ ml. 


(14) 


Consider any nonempty collection of columns of A denoted by C . Take Ci = C and C2 = 0. Then the sum 
of the columns in Ci minus the sum of the columns in C 2 will be a vector with entries 0,1 and —1, because 
in A there cannot exist more than one row corresponding to an arc between two nodes of the constraint 
graph and each such row has exactly two nonzero entries, a 4-1 and a —1. So, by Theorem[l]the polyhedron 
{x G RI^I : I ^ Ax A w, 0 ^ a: ^ ml} has only integral vertices and optimizing the linear objective in 
problem (14) over this polyhedron will result in an integral solution. □ 
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After solving the linear programming problem (13), we obtain an integral timetable, which we will call 
the energy minimizing timetable (EMT). We denote the optimal decision vector of this timetable by x in 
the constraint graph notation and in the original notation. 


5 Second optimization model 

In this section we modify the trip time constraints such that the total energy consumption of the final 
timetable is kept at the same minimum as the EMT. Then, we describe our optimization strategy aimed to 
maximize the utilization of regenerative energy of braking trains, and we present the second optimization 
model. 


5.1 Keeping the total energy consumption at minimum 

In any feasible timetable, if the trip times are kept to be the same as the ones obtained from the EMT, 
then the energy optimal speed profiles for all trains will be the same. As a result, the energy consumption 
associated with that timetable will remain at the same minimum as found in the EMT. So, in the second 
optimization problem, instead of using the trip time constraint, for every trip we fix the trip time to the 
value in the EMT, i.e., 


yteT, 

y{i,j)eA*, a]-d\=a]-d\, 

(15) 

V(b j) e <p, 

V(t, t') e Bij, a*' -d\= a*' - d\. 

(16) 


For all other constraints, bounds are allowed to vary as described by Equations (§ §. As a consequence of 
fixing all trip times, the power graph of every trip made by any train becomes known to us, since it depends 
on the corresponding optimal speed profile calculated in real-time by existing software US], in page 285]. 

5.2 Maximizing the utilization of regenerative energy of braking trains 

In this subsection we describe our strategy to maximize the utilization of the regenerative energy produced 
by the braking trains. Strategies based on transfer of regenerative braking energy back to the electrical grid 
requires specialized technology such as reversible electrical substations [55j. A strategy based on storing is 
not feasible with present technology, because storage options such as super-capacitors, fly-wheels, e.t.c., have 
drastic discharge rates besides being too expensive [29l page 66], [30l page 92]. A better strategy that can 
be used with existing technology [31] is to transfer the regenerative energy of a braking train to a nearby 
and simultaneously accelerating train, if both of them operate under the same electrical substation. We call 
such pairs of trains suitable train pairs. So our objective is to maximize the total overlapped area between 
the graphs of power consumption and regeneration of all suitable train pairs. To model this mathematically, 
we are faced with the following tasks: i) define suitable train pairs, ii) provide a tractable description of the 
overlapped area between power graphs of such a pair. We describe them as follows. 

5.2.1 Defining suitable train pairs 

We consider platform pairs who are opposite to each other and are powered by the same electrical substations. 
Thus, the transmission loss in transferring electrical energy between them is negligible. In our work we The 
set containing all such platform pairs is denoted by fl. Consider any such platform pair {i,j) S and let 
7i C 7” be the set of all trains which arrive at, dwell and then depart from platform i. Suppose, t G %. 
Now, we are interested in finding another train t on platform j, i.e., t € Tj, which along with t would form 
a suitable pair for the transfer of regenerative braking energy. To achieve this, we use the EMT. Among all 
trains going through platform j, the one which is temporally closest to t in the energy-minimizing timetable 
is be the best candidate to form a pair with t. The temporal proximity can be of two types with respect to 
t, which results in the following definitions. 
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Definition 1. Consider any {i,j) G fi. For every train t G Ti, the train t G Tj is called the temporally 
closest train to the right of t if 


t = 


argmin 




^<r} 



(17) 


where r is an empirical parameter determined by the timetable designer and is much smaller than the time 
horizon of the entire timetable. 

Definition 2. Consider any {i,j) G fl. For every train t G Ti, the train t G 7} is called the temporally 
closest train to the left of t if 




t'e{xeTj-.o<- 


+ d^ 


-^<rl 



(18) 
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Definition 3. Consider any {i,j) G fi. For every train t G %, the train t G Tj is called the temporally 
closest train to t if 

t = argmin 

If both t and t are temporally equidistant from t, we pick one of them arbitrarily. 

Any synchronization process between a suitable train pair (SPSTP) can be described by specifying the 
corresponding i, j, t and i by using the definitions above. We construct a set of all SPSTPs, which we denote 
by E. Each element of this set is a tuple of the form {i,j, t, t). Because t is unique for any t in each element 

of £, we can partition £ into two sets denoted by £ and £ containing elements of the form t ) and 

t) respectively. 

Remark. Now we compare our design choice to that of El, which considers two additional scenarios: two 
trains 1) on two consecutive platforms on the same track, or 2) on the same platform on same track. In the 
first case, there will of course be a significant distance between the two trains, and in the second case, we 
need to maintain at least the safety distance between the two trains irrespective of how small the headway 
is. First, to transfer the regenerative energy to the accelerating trains, we would need either supercapacitors, 
or fly-wheels, both of which have drastic discharge rates [29l page 66], [30l page 92]. Additionally, if we 
consider resistive transmission loss, in practice there would be very little or no energy transfer due to the 
power loss along the line. Note that in m this issue would not arise, as they have assumed no transmission 
loss. 



5.2.2 Description of the overlapped area between power graphs 

The power graph during accelerating and braking is highly nonlinear in nature with no analytic form, as 
shown in Figure So, maximizing the exact overlapped area will lead to an intractable optimization 
problem. However, the existence of dominant peaks with sharp falls allows us to apply a robust lumping 
method such as j heuristic |32l page 33-34] to approximate the power graphs as rectangles. The ^ heuristic 
is applied as follows (see Figure]^. The height of the rectangle is the maximum power, and the width is 
the interval with extreme points corresponding to power dropped at 1/e of the maximum. For the sharp 
drop from the peak, such rectangles are very robust approximations to the original power graph containing 
the most concentrated part of the energy, e.g.,, if the drop were exponential, then the energy contained by 
the rectangle would have been exactly equal to that of the original curve (32 page 33-34]. After converting 
both the power graphs to rectangles, maximizing the overlapped area under those rectangles is equivalent to 
aligning the midpoint of the width of the rectangles; we call such a midpoint regenerative or consumptive 
alignment point. These alignment points act as virtual peaks of the approximated power graphs. As shown 
in Figure]^ for a train t in its braking phase prior to its arrival at platform *, the relative distance of a* 
from the regenerative alignment point is denoted by V*, while during acceleration the relative distance of 
the consumptive alignment point from d* is denoted by A*. Note that both relative distances are known 
parameters for the current optimization problem. 


5.3 Second optimization model 


Consider an element t) G £. To ensure the transfer of maximum possible regenerative energy from 

the braking train t to the accelerating train t, we aim to align both their alignment points such that 
dl+ A*= Uj — Vj, or keep them as close as possible otherwise. Similarly, for any t) G £, our 

objective is d- 


Aj* = a* — V* , or as close as possible. Let a decision vector y be defined as 


y = 


(M 




^ t)e£ ^ ^ ^ t )e£ 


( 20 ) 


15 






Then our goal comprises of two parts: 1) maximize the number of zero components of y which corresponds 
to minimizing card( 2 /), and 2 ) keep the nonzero components as close to zero as possible which corresponds 
to minimizing the norm of y, ||j/||i. Combining these two we can write the exact optimization problem as 
follows: 

minimize card( 2 /) + 7 ||j/||i 
subject to 

Equations^ - ([^, ([^, ( [^ , 

0 < a- < TO, 0 < ^ < m, Vi G Af*, Vt G T, 


( 21 ) 


where 7 is a positive weight, and decision variables are a, d and y. The objective function is nonconvex as 
shown next. Take the convex combination of the vectors 2 ei /7 0 with convex coefhcients 1/2. Then, 


card ( — I + 7 

7 




2ei 

7 


i (card (0) +7||0|li) = 1.5, 


and thus violates definition of a convex function. As a result, problem (21) is a nonconvex problem. Note 


that if we remove the cardinality part from the objective, then it reduces to a convex optimization problem 
because the constraints are affine and the objective is the ii norm of an affine transformation of the decision 
variables [33l pages 72, 79, 136-137]. Such problems are often called convex-cardinality problem and are of 
AfT^-hard computational complexity in general |34j . An effective yet tractable numerical scheme to achieve a 
low-cardinality solution in a convex-cardinality problem is the £i norm heuristic, where card( 2 /) is replaced 
by II2/11 1 j thus converting problem ( 21 ) into a convex optimization problem. This is described by problem ( 22 ) 


below. The ii norm heuristic is supported by extensive numerical evidence with successful applications to 
many fields, e.g.,, robust estimation in statistics, support vector machine in machine learning, total variation 
reconstruction in signal processing, compressed sensing etc. In the next section we show that in our problem 
too, the £i norm heuristic produces excellent results. Intuitively, the ii norm heuristic works well, because it 
encourages sparsity in its arguments by incentivizing exact alignment between regenerative alignment points 
with the associated consumptive ones |3S1 pages 300-301]. We provide a theoretical justification for the use 
of £1 norm in our case as follows. 

Proposition 2. The convex optimization problem described by 


minimize 
subject to 
Equations([3l) — ([^, ((l5|, ([l 6 |, ([2^ 

n ^ ^ n ^ ^ ^ ..nr, VIn r Kf 


( 22 ) 


0 < a- < TO, 0 < ^ < TO, Vi G Af*, Vt G T, 


is the best convex approximation of the nonconvex problem (21) from below. 


Proof. Both problems (22) and (21 1 have the same constraint set, so we need to focus on the objective only. 


The best convex approximation of a nonconvex function / : C — R (where C is any set) from below is 
given by its convex envelope env / on C. The function env / is the largest convex function that is an under 
estimator of / on C, i.e., 

env / = sup{/ : C — R I / is convex and / < /}, 

where sup stands for the supremum, , the least upper bound of the set. The definition implies, epienv/ = 
conv epi /. 


From Equation (20) we see that y is an affine transformation of a and d, and from the last constraints of 
problem (21) we see that both a and d are upper bounded by to, i.e., ||a||oo < m and ||d||oo < w- So there 
exists a positive number P such that ||2/||oo < P- As the domain of y is bounded in an i^o ball with radius 
P, envcard(?/) = 4||2/||i [Ml page 321]. So, the best convex approximation of the objective from below is 
;pll2/lli+7ll2/lli = (p+7)ll2/lli- As the coefficient (^ -I- 7 ) is a constant for a particular optimization problem, 
it can be omitted, and thus we arrive at the claim. □ 
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Using the epigraph approach [531 pages 143-144], we can transform the convex problem (22) into a linear 
program as follows. For each t) G £ and each t) G £, we introduce new decision variables 

and respectively, such that 0“ > |c?--|- A- —aj + | and 0“ > \dj + —a- -I- V-|. Then, the 

convex optimization problem can be converted into the following linear problem: 


minimize ^ ^ 0*/ 

subject to 

dU > dl+ A* -o/ -f V *, 


Qtt 


3 3 _ 

't At , „t V7 t 


,3 > -d\- A* +a/ - V] 


> d* + A * —a- -I- V', 


Qt t 

> -d* - a/ +al - V‘ 




t)G £^ 
t) Gj. 
V(i,j,t, t)G£^ 
t) G £ 


Equations (|^ — ([s]), (15), (16), 

0 < a* < m, 0 < (^< m yt gT, Vj S Af*, 


(23) 


where the decision variables are a■, d■, 6*) * and d 


6 Limitations of the optimization model 

In this section, we discuss the limitations of our model as follows. 

• We have assumed that the amount by which the trip time is allowed to vary is on the order of seconds 
(Assumption [^. Though this is true for most of the subway systems, there are exceptions where 
this assumption may not hold. For example, when a trip between two cities is considered (especially 
involving different countries), the trip is on the order of hours with the acceptable trip time bound 
often being on the order of 5-10 minutes and even more in some cases. In such a scenario, an affine 
approximation of the energy with respect to trip time would not be very efficient any more, and our 
model would not be suitable for such a case. 

• In Section [^ we have have applied two different heuristics to arrive at a convex optimization problem. 
At first we have used ^ heuristic to come up with a tractable description of the overlapped area 
between power graphs, and then we have used the ii norm heuristic to approximate a nonconvex 
objective function with its convex envelope. So, it is quite likely that the timetable obtained by 
solving the convex optimization problem (Problem |23| ) would have a worse objective value compared 
to the original intractable optimization problem. For this reason our model is energy-efficient, but not 
necessarily energy-optimal. 

• The model does not directly address the case of significant delay. However, we have considered two 
indirect ways of dealing with it in practice. 

— In any automatic train supervision system, which has the responsibility of implementing the 
timetable, dwell and velocity regulation are performed to maintain trains on their proper time. If 
there is a deviation from the optimal timetable because of some delay, the ATS performs regulation 
to move delayed trains back to the planned optimal timetable timetable. Thus the system will 
typically return to a normal state in less than half an hour after a delay of one minute. 

— Another way is incorporating the change in the system (due to the delay) as an input data and 
solving a new bnt shorter optimization problem with a time horizon of 1-2 hours which can be 
done in real time using our model. While the shorter model is being implemented we solve the 
larger optimization problem spanning the rest of the day. 
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Figure 7: Railway network considered for numerical study 


7 Numerical study 

In this section we apply our model to different problem instances spanning full service period of one day to 
service PES2-SFM2 of line 8 of Shanghai Metro network. Shanghai Metro is the world’s largest rapid transit 
system by route length, second largest by number of stations(after Beijing), and third oldest rapid transit 
system in mainland China. Line 8, opened on December of 2007, is one of the 14 lines of the Shanghai Metro 
Network. It passes by some of Shanghai’s densest neighborhoods, and has a daily ridership of approximately 
1.08 million (2014 peak). This line is 37.4 km long with 28 stations in operation [3^. The service PES2- 
SFM2 of line 8 of the Shanghai Metro network is shown in Figure There are two lines in this network: 
Line 1 and Line 2. There are fourteen stations in the network denoted by all capitalized words in the figure. 
Each station has two platforms each on different train lines, e.g.,, LXM is station that has two opposite 
platforms: LXMl and LXM2 on Line 1 and Line 2 respectively. The platforms are denoted by rectangles. 
The platforms indicated by PES2 and SFM2 are the turn-around points on Line 2, with the crossing-overs 
being PES2-GRW1 and LHS1-SFM2. The number of trains, headway times, speed of the involved trains, 
the grades of the tracks and nature of the energy profile of the associated acceleration and braking phases 
are different for different instances. 

Now we provide some relevant information regarding the railway network in consideration. The line is 37 
km long. The average distance between two stations in this network is 1.4 km, with the minimum distance 
being 738 m (between YHR and ZJD) and maximum distance being 2.6 km (between PJT and LHR). The 
slope of the track is in [—2.00453°, 2.00453°]. The maximum allowable acceleration of a train is 1.04 m/s 
at accelerating phase. The maximum allowable deceleration rate of a train during coasting phase is -0.2 
m/s^. The maximum allowable deceleration of a train during braking phase is -0.8 m/s^. The conversion 
factor from electricity to kinetic energy is 0.9, and the conversion factor from kinetic to regenerative braking 
energy is 0.76. The transmission loss factor of regenerative electricity is 0.1. The mass of the train is in 
[229370,361520] kg, with the average mass being 295445 kg. 

The data on speed limit is described by Table in Appendix. The speed limit data are based on grade 
and curvature of the tracks, and operational constraints present in the system. For the railway network, the 
tracks have piece-wise speed limits for trains, i.e., each track (except CSR2-YHR2) is divided into multiple 
segments, where each segment has a constant speed limit. The track CSR2-YHR2 has only one segment 
(itself) with a speed limit of 60 km/h. Table is provided to us by Thales Canada Inc. 

The numerical study was executed on a Intel Core 17-46400 CPU with 8 GB RAM running Windows 8.1 
Pro operating system. For modelling the problem, we have used JuMP - an open source algebraic modelling 
language embedded in programming language Julia [36] . Within our JuMP code we have called academic 
version of Gurobi Optimizer 6.0 as the solver. We have implemented an interior point algorithm because of 
the underlying sparsity in the data structure. As mentioned before, a measure of the quality of affine fittings 
using least-squares approach is given by the coefficient of determination^ which can vary between 0 to 1, 
with 0 being the worst and 1 being the best EH page 518]. In our numerical study, the average coefficient of 
determination of the affine fittings for energy versus trip times over all different trips and all trains is found 
to be 0.9483 with a standard deviation of 0.05, which justifies our approach. 

The duration of the timetables is eighteen hours which is the full service period of the railway network. 
We have considered eleven different instances with varying average headway times and number of trains. 
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The number of trains increases as the average headway time decreases, where the relation between them can 
be determined from Equation The results of the numerical study are shown in Table 

We can see that, in all of the cases our model has found the optimal timetables very quickly, the largest 
runtime being 12.58s. To the best of our knowledge, this model is the only one to calculate energy-efficient 
railway timetable spanning an entire day, the next largest being 6 hours only [5] with a much larger com¬ 
putation time for smaller sized problems. After we get the final timetable, we calculate the total effective 
energy consumption by all trains involved in SPSTPs and compare it with the original timetables. The 
effective energy consumption of a train during a trip is defined as the difference between the total energy 
required to make a trip and the amount of energy that is being supplied by a braking train during syn¬ 
chronization process. So, the effective energy consumption is the energy that will be consumed from the 
electrical substations. 

Now we briefly point out how the total energy consumption is calculated theoretically. Consider the trip 
of any train t of mass mt from one platform to the next one. Assume the trip time to make the trip is T. At 
time instance r, the acceleration, speed, position, and resistive acceleration is denoted by Ut{T),Vt{T), Stir) 
and —r (vtir))^ respectively. The resistive acceleration is given by the Davis formula [22] 

^ (^t(T)) = ao + aiVtir) a2V^{T), 

where Oq, Oi, 02 are known positive numbers. At time instance r, the net acceleration of the train is equal to 
M((r) — r (^((t)), where ^((t) bounded by magnitude constraints which are non-increasing and depends on 
the speed of the train. Only the positive acceleration force consumes energy and negative acceleration during 
braking produces regenerative energy. Define, u^(t) = max{0, Ut(r)}, and uf {t) = max{0, —it(T)} which 
represents the positive and negative part of Ut^r), respectively. The negative acceleration associated with 
regenerative braking is denoted by Ut”regen('^)' the total energy consumed by the train while making this 

trip is given by the integral: mtuf {T)st{T)dT and the total regenerative energy produced by the train is 

given by mtuf^^^^^(T)st(T)dT. 

Consider a suitable train pair, say trains t and t associated with the SPSTP {i,j,t,t), with train t accel¬ 
erating and train t braking. After we have calculated the timetable, it will provide us with the duration time 
during which they are synchronized. Let us denote the times for the beginning and end of the synchroniza¬ 
tion process by Ti and T 2 . For any x G R, denote x+ = max{x,0}. Then the effective energy consumption 
of associated with this train pair is denoted by, 

(mtut (T)st(r) - mtUf^^^^^jT)si{T)^ dr, 

The physical interpretation for the above integral is as follows. If for some reason we have more regenerative 
energy can be provided than needed by accelerating train (though the likelihood of the occurrence of this 
case is very low in practice), then the extra energy is burned via resistive braking. Similarly we can define 
the effective energy consumption with train t accelerating and train t braking as follows: 

/-T2 

/ {miU^{T)si{T)-mtuf,^g^^{T)st{T)) dr. 

JTi 



So, the effective energy consumption effective energy consumption associated with the SPSTP {i,j,t,i) is 
denoted by: 


E. 


effective 



mtuf {T)st{T) — rnfU~ train t is accelerating and train i is braking 


{t)si{t) — rntUf j.ggg^{T)st{T)'^ dr, if train t is braking and train i is accelerating 


In practice, these integrations are performed numerically for which robust and fast packages exist e.g., 
in our case. The total effective energy consumption over all SPSTPs is given by: 


E, 


effective 


■^effective’ 
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Table 2: Results of the numerical study performed to line 8 of Shanghai Metro network 


Number 

of 

trains 

Number 

of con¬ 
straints 

step 1 

Number 

of vari¬ 
ables 

step 1 

step 1 

CPU 

time (s) 

Number 

of con¬ 
straints 

step 2 

Number 

of Vari¬ 
able 

step 2 

step 2 

CPU 

Time 

(s) 

Initial effec¬ 
tive energy 

consumption 

(kWh) 

Final effec¬ 
tive energy 

consumption 

(kWh) 

Reduction in 

effective energy 

consumption 

1000 

91998 

30060 

3.24 

116558 

34871 

6.03 

250951.3 

201658.7 

19.64 % 

1032 

94944 

31022 

3.03 

120394 

36038 

5.45 

261994.5 

208558.1 

20.40 % 

1066 

98074 

32044 

3.96 

124494 

37290 

5.62 

272486.7 

215896.7 

20.77 % 

1100 

101204 

33066 

3.47 

129284 

38887 

5.39 

288677.5 

229091.8 

20.64 % 

1132 

104150 

34028 

3.15 

133354 

40171 

6.69 

308924.4 

243672 

21.12 % 

1166 

107280 

35050 

2.84 

137322 

41357 

6.67 

322612.7 

256288.7 

20.56 % 

1198 

110226 

36012 

2.96 

141322 

42606 

7.16 

329388.2 

262205.6 

20.40% 

1232 

113356 

37034 

4.04 

145756 

44025 

7.61 

354050.2 

277536.7 

21.61% 

1266 

116486 

38056 

3.84 

149868 

45283 

8.74 

368901.4 

297815 

19.27 % 

1298 

119432 

39018 

3.93 

153480 

46338 

7.62 

366488.4 

293068.8 

20.03 % 

1332 

122562 

40040 

4.22 

157752 

47676 

8.02 

379700.8 

300910.1 

20.75 % 


The original timetables, which we compare the final timetables with, are provided by Thales Canada Inc. 
It should be noted that, the number of trains T is fixed for each of the instances. The energy calculation is 
done using SPSIM, which is a proprietary software owned by Thales Canada Inc m, and Cubature, which 
is an open-source Julia package written by Steven G. Johnson that uses an adaptive algorithm for the 
approximate calculation of multiple integrals [35] ■ SPSIM calculates the power versus time graphs of all the 
active trains for the original and optimal timetables. Cubature is used to calculate the effective area under 
the power versus time graphs to determine 1) the total energy required by the trains during the trips, 2) 
the total transferred regenerative energy during the SPSTPs, and 3) the effective energy consumption as 
the difference of the first two quantities. The effective energy consumption of the optimal timetables in 
comparison with the original ones is reduced quite significantly - even in the worst case, the reduction in 
effective energy consumption is 19.27%, with the best case corresponding to 21.61%. 


8 Conclusion 

In this paper we have proposed a novel two-step linear optimization model to calculate an energy-efficient 
timetable in modern metro railway networks. The objective is to minimize the total electrical energy con¬ 
sumption of all trains and to maximize the utilization of regenerative energy produced by braking trains. In 
contrast to other existing models, this model is computationally the most tractable one. We have applied 
our optimization model to eleven different instances of service PES2-SFM2 of line 8 of Shanghai Metro 
network. All instances span the full service period of one day (18 hours) with thousands of active trains. 
For all instances our model has found optimal timetables in less than I3s with significant reductions in the 
effective energy consumption. Code based on our optimization model has been integrated with the industrial 
timetable compiler of Thales Inc. 
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Appendix 


Table 3: Speed limit for line 8 of Shanghai Metro network 


Origin-Destination 

Start (m) 

End (m) 

Speed limit (km/h) 

CSRl-YSSl 

0.0 

143.5 

60 


143.5 

1004.6 

70 


1004.6 

1138.2 

60 

CSR2-YHR2 

0.0 

910.0 

60 

GRWl-LXMl 

0.0 

153.3 

60 


153.3 

870.1 

70 


870.1 

1006.9 

60 

GRW2-PES2 

0.0 

173.1 

60 


173.1 

636.4 

70 


636.4 

769.5 

60 

JYRl-LZVl 

0.0 

1366.7 

60 


1366.7 

2220.6 

65 


2220.6 

2357.3 

60 

JYR2-YSS2 

0.0 

143.4 

60 


143.4 

1388.9 

70 


1388.9 

1522.3 

60 

JYSl-LHSl 

0.0 

140.0 

60 


140.0 

829.2 

75 


829.2 

1202.3 

60 

JYS2-PJT2 

0.0 

140.0 

60 


140.0 

371.1 

70 


371.1 

1081.9 

75 


1081.9 

1249.8 

70 


1249.8 

1386.2 

60 

LHRl-PJTl 

0.0 

140.1 

60 


140.1 

766.2 

70 


766.2 

1623.4 

75 


1623.4 

1805.9 

70 


1805.9 

2374.4 

75 


2374.4 

2487.8 

70 


2487.8 

2624.3 

60 

LHR2-LZV2 

0.0 

139.8 

60 


139.8 

2457.3 

70 


2457.3 

2594.1 

60 

LHSl-SFMl 

0.0 

140.0 

60 


140.0 

1220.1 

70 


1220.1 

1357.4 

60 

LHS2-JYS2 

0.0 

186.7 

60 


186.7 

853.8 

75 


853.8 

1064.4 

70 


1064.4 

1200.8 

60 

LJBl-SXZl 

0.0 

140.1 

60 
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YHRl-CSRl 

0.0 

143.4 

60 


143.4 

773.7 

70 


773.7 

910.3 

60 

YHR2-ZJD2 

0.0 

140.0 

60 


140.0 

601.3 

70 


601.3 

738.0 

60 

YSSl-JYRl 

0.0 

140.2 

60 


140.2 

664.7 

70 


664.7 

987.7 

75 


987.7 

1389.2 

70 


1389.2 

1525.9 

60 

YSS2-CSR2 

0.0 

430.4 

60 


430.4 

1014.8 

70 


1014.8 

1151.6 

54 

ZJDl-YHRl 

0.0 

140.1 

60 


140.1 

605.9 

70 


605.9 

742.5 

60 

ZJD2-SXZ2 

0.0 

353.8 

60 


353.8 

1393.7 

70 


1393.7 

1910.0 

65 


1910.0 

2043.3 

60 
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